#include<opencv2/opencv.hpp>
#include<string>
#include<vector>
#include<Eigen/Sparse>
#include<cmath>
#include<iostream>
#include<iomanip>
using std::string;
namespace ei=Eigen;
#define square(x) ((x)*(x))

constexpr int W=100,H=100;
#undef ddot
auto sampling(const string& path)->cv::Mat{
    cv::Mat img = cv::imread(path,cv::IMREAD_COLOR);
    cv::GaussianBlur(img,img,cv::Size(11,11),11);
    cv::resize(img,img,cv::Size(100,100));
    return img;
}

auto Laplacian_mat(const cv::Mat& inimg)->ei::SparseMatrix<double>{
    int count=0,used=0;
    cv::Mat img=inimg.clone();
    img.convertTo(img,CV_64FC3);
    img/=255.0;
    std::vector<ei::Triplet<double> > triplets;
    ei::VectorXd D(10000);
    for (int i = 0; i < 100; i++){
        for (int j = 0; j<100; j++){
            double d=0;
            for (int k = W-2; k <= W+2; k++){
                for (int l = H-2; l <= H+2; l++){
                    int kk=(k>=W?k-W:k);
                    int ll=(l>=H?l-H:l);
                    cv::Vec3d diff=img.at<cv::Vec3d>(i,j)-img.at<cv::Vec3d>(kk,ll);
                    double dist=diff.dot(diff)+pow((i-kk)/100.0,2)+pow((j-ll)/100.0,2);
                    if(!used){
                        std::cout<<dist<<'\n';
                        used=true;
                    }
                    dist=exp(-4 * dist);
                    triplets.emplace_back(100*i+j,100*kk+ll,dist);
                    ++count;
                    d+=dist;
                }
            }
            if (d>0){
                D[100*i+j]=pow(d,-0.5);
            } else {
                D[100*i+j]=0;
            }
            
            
        }
    }

    std::cout<<"Total edges:"<<count<<'\n';
    std::cout<<"average degrees:"<<((double)count/(double)(W*H))<<'\n';
    ei::SparseMatrix<double> G(W*H,W*H);
    G.setFromTriplets(triplets.begin(),triplets.end());
    G.makeCompressed();
    ei::SparseMatrix<double> diag(D.asDiagonal());
    ei::SparseMatrix<double> I(W*H,W*H);
    I.setIdentity();
    ei::SparseMatrix<double> L=I+diag*G*diag;
    L.makeCompressed();
    std::cout<<"coeff=" << L.coeff(9898, 0) << '\n';

    return L;
}

auto Getting_eigen_vec(const ei::SparseMatrix<double>& L)->ei::VectorXd{
    constexpr int k=500;
    ei::VectorXd vec=ei::VectorXd::Random(L.cols())+ei::VectorXd::Ones(L.cols());
    vec/=2;

    for(int i=0;i<k;i++){
        vec=L*vec;
    }

    vec.normalize();

    return vec;
}

auto kmeans(ei::VectorXd vec)->cv::Mat{
    vec*=10000.0;
    cv::Mat img(cv::Size(W,H),CV_8U);
    img=0;
    double m1=vec[0],m2=vec[W/2 * 100 + H/2];
    constexpr int k=10;
    for (int _ = 0; _< k; _++){
        double m1_nxt=0,m2_nxt=0;
        int count1=0,count2=0;
        for (int i = 0; i < W; i++){
            for(int j = 0; j < H; j++){
                if(square(vec[100*i+j]-m1) <= square(vec[100*i+j]-m2)){
                    img.at<uchar>(i,j)=0;
                    m1_nxt+=vec[100*i+j];
                    count1+=1;
                }else{
                    img.at<uchar>(i,j)=255;
                    m2_nxt+=vec[100*i+j];
                    count2+=1;
                }
            }
            if(count1==0) count1=1;
            if(count2==0) count2=1;
            m1=m1_nxt/count1;
            m2=m2_nxt/count2;
        }
        
    }
    
    std::cout<<vec[0]<<' '<<vec[5050]<<'\n';
    return img;
}

int main(int argc, char** argv){
    std::cout<<std::fixed<<std::setprecision(19);
    string path;
    if(argc==1){
        path="E:/repos/imgspec/star.jpg";
    } else {
        path=argv[1];
    }

    cv::Mat img = sampling(path);

    ei::SparseMatrix<double> L=Laplacian_mat(img);
    
    ei::VectorXd x=Getting_eigen_vec(L);
    printf_s("Lcol=%lld\n", L.cols());
    cv::Mat img_new=kmeans(x);

    cv::namedWindow("split picture", cv::WINDOW_NORMAL);
    cv::namedWindow("Downsampled", cv::WINDOW_NORMAL);
    cv::resizeWindow("split picture", 400, 400);
    cv::resizeWindow("Downsampled", 400, 400);

    cv::imshow("Downsampled", img);
    cv::imshow("split picture", img_new);
    cv::waitKey(0);
}



